function [out_out,DAT_Trace,FractionatingMode_Out,bulkD,PLAG_comp,CHANGEINDICIES] = YKG_wTRACE_Magnetite(inputfrommodel,PINC,trace,Dmatrix)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This mfile is used to in conjunction with Kinzler and Grove 1992.% it is an adaptation of Yang et al. 1996 from FORTRAN to MATLAB by PM % Gregg 2007.  It calculates fractional crystallization paths given% an initial mantle composition.  This version of the code iterates% to look at compositions from 4 kbar to 0.001 kbar%% INPUT:    input - mantle composition.  must be 12 component and %               in the following order: %               SIO2, TIO2, AL2O3, FE2O3, FEO, MGO, CAO, NA2O, K2O,P2O5,%               CR2O3, MNO%           outfilenm - this is a character string for the output file%               which contains the oxide wt of the residual melt for each%               fractional xlization step% Outputs can vary as needed.% OUTPUT:   In the current version of the code outputs are all written to%           output files.  These files can vary based on what information%           the user wants outputed.% %% T. Gregg 2007%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%input = [46.58 0.95 17.66 0 9.31 11.90 11.43 1.81 0.37 0 0 0];%input = [48.62 0.94 18.31 0 7.92 10.09  9.93 3.82 0.37 0 0 0];%input = [49.37 0.89 17.49 0 7.98 10.28 10.46 3.31 0.23 0 0 0];% inputfrommodel% [SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,CaO,K2O,Na2O] wt%;inputfrommodel(isnan(inputfrommodel))=0; inputfrommodel = inputfrommodel./nansum(inputfrommodel)*100;input = inputfrommodel; %input = [inputfrommodel(1:3) 0 inputfrommodel(5:7) inputfrommodel(9) inputfrommodel(8) 0 inputfrommodel(4) 0];out = [];% KDOL = 0.3;% KDCPX = 0.23;% KDOL = 0.3;% KDCPX =.94*KDOL;%0.2726;%KDOL = 0.23;KDOL = 0.3;KDCPX = 0.23;%PINC = [1];%Main output file%file = ['runs-05.11.09/R2Du5Nu41e23_refinedT1325WH_YangOut_04kbar'];%Find Fe8 and Na8file2 = ['runs-05.11.09/Na8Fe8_out'];    bulkD = []; kk = 0;for j = 1:1%fid6 = fopen([file2(j,:),'.m'],'wt');for i = 1:length(PINC)    i;    Mg8cnt = 1;    kk = kk+1;    DAT = zeros(5,12);    DAT(1,1:12) = input(j,:);    DAT_Trace(1,:) = trace;    FractionatingMode_Out(1,:) =[0 0 0 0 0 0 0];         PLAG_comp = zeros(1,12);    OL_comp= zeros(1,12);    %fid8 = fopen([file(kk,:),'.m'],'wt');    P = PINC(i);    %fprintf(fid8,['Out= [']);    %fprintf(fid6,'P(');fprintf(fid6,'%2.0f',i);fprintf(fid6,') = ');    %fprintf(fid6,'%8.4f',P);fprintf(fid6,';');    CHANGEINDICIES = [];     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    %  CRYSTAL FRACTIONATION OR ADDITION    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    CHECK1 = 0;    %CHECK1 = input(['\nENTER 0: FOR CRYSTAL FRACTIONATION \n',...    %     '      1: FOR OL-PL-AUG ADDITION. \n',...    %     '\n NOTE: THIS PROGRAM DOES NOT PERFORM OL OR OL-PL ADDITION. \n',...    %     'HUAI-JEN YANG HAS ANOTHER PROGRAM FOR DOING THESE. JUST ASK.\n']);    % (CHECK1.EQ.1) GOTO 500    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    %  EVALUATING OLIV OR PLAG SATURATED    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    [QTZ,PLAG,OLIV,CPX,QQTZ,QOLIV,QCPX,PPLAG,POLIV,PCPX] = Tormey(DAT(1,:));    [ROQTZ,ROOLIV,ROPLAG,ROCPX,PROPLAG,PROOLIV,PROCPX,QROQTZ,...        QROOLIV,QROCPX] = OPALM_TRACE(DAT,P);    [XOPALM,YOPALM] = PLOT_TERN_TRACE(PROOLIV,PROCPX);    [XMODE,YMODE] = PLOT_TERN_TRACE(POLIV,PCPX);    DIFF1=XOPALM-XMODE;         test = 1;    if DIFF1 > 0        %%% 100        % Olivine fractionation        while (DIFF1 > 0)             [DAT(2,1),DAT(2,5),DAT(2,6)] = OLCOM_TRACE(DAT,KDOL);                     for J=1:12                DAT(1,J)=(DAT(1,J)-DAT(2,J)*0.01)/0.99;            end            FractionatingMode = [0 0 1 0 0 0 0];                 D = FractionatingMode*Dmatrix';              bulkD(:,test+1) = [D];            %D = Dmatrix*FractionatingMode';                    DAT_Trace(test+1,:) = (DAT_Trace(test,:) - 0.01.*(DAT_Trace(test,:).*D))./0.99;             FractionatingMode_Out(test+1,:) =FractionatingMode;                         out = [out;DAT(1,1:12)];            %fprintf(fid8,'%8.4f',DAT(1,1:12));            %fprintf(fid8,' \n');            [QTZ,PLAG,OLIV,CPX,QQTZ,QOLIV,QCPX,PPLAG,POLIV,PCPX] = Tormey(DAT(1,:));            [ROQTZ,ROOLIV,ROPLAG,ROCPX,PROPLAG,PROOLIV,PROCPX,QROQTZ,...                QROOLIV,QROCPX] = OPALM_TRACE(DAT,P);            [XOPALM,YOPALM] = PLOT_TERN_TRACE(PROOLIV,PROCPX);            [XMODE,YMODE] = PLOT_TERN_TRACE(POLIV,PCPX);            DIFF1=XOPALM-XMODE;            test = test+1;            end                %SMB: is this because it does an extra step?  should the previous        %step just be deleated? %         for J=1:12%             DAT(1,J)=(DAT(1,J)+DAT(2,J)*0.01)/1.01;%         end%         out = [out;DAT(1,1:12)];        %test        %fprintf(fid8,'\n');        %fprintf(fid8,'%8.4f',DAT(1,1:12));    else        %%% 21        % Plagioclase fractionation               while(DIFF1 < 0)                        [DAT(3,1),DAT(3,3),DAT(3,7),DAT(3,8)] = OPAMPL_TRACE(DAT,P);                        for J=1:12                DAT(1,J)=(DAT(1,J)-DAT(3,J)*0.01)/0.99;            end            out = [out;DAT(1,1:12)];                                    FractionatingMode = [0 0 0 1 0 0 0];                        D = FractionatingMode*Dmatrix';            bulkD(:,test+1) = [D];            %D = Dmatrix*FractionatingMode';                     DAT_Trace(test+1,:) = (DAT_Trace(test,:) - 0.01.*(DAT_Trace(test,:).*D))./0.99;            FractionatingMode_Out(test+1,:) =FractionatingMode;                                     %fprintf(fid8,'\n');            %fprintf(fid8,'%8.4f',DAT(1,1:12));            [QTZ,PLAG,OLIV,CPX,QQTZ,QOLIV,QCPX,PPLAG,POLIV,PCPX] = Tormey(DAT(1,:));            [ROQTZ,ROOLIV,ROPLAG,ROCPX,PROPLAG,PROOLIV,PROCPX,QROQTZ,...                QROOLIV,QROCPX] = OPALM_TRACE(DAT,P);            [XOPALM,YOPALM] = PLOT_TERN_TRACE(PROOLIV,PROCPX);            [XMODE,YMODE] = PLOT_TERN_TRACE(POLIV,PCPX);            DIFF1=XOPALM-XMODE;                        test = test+1;        end               %%         for J=1:12%             DAT(1,J)=(DAT(1,J)+DAT(3,J)*0.01)/1.01;%         end    end    %% olivine and plag fractionation    A=100;    %%% 200    B = 1;    %test = 0;    OP_first = 0;     while B > 0;                          [DAT(2,1),DAT(2,5),DAT(2,6)] = OLCOM_TRACE(DAT,KDOL);	            [DAT(3,1),DAT(3,3),DAT(3,7),DAT(3,8)] = OPAMPL_TRACE(DAT,P);                        [ROQTZ,ROOLIV,ROPLAG,ROCPX,PROPLAG,PROOLIV,PROCPX,QROQTZ,...            QROOLIV,QROCPX] = OPALM_TRACE(DAT,P);                [XOPALM,YOPALM] = PLOT_TERN_TRACE(PROOLIV,PROCPX);                PLPROP=XOPALM/2*1.732;        OLPROP=1-PLPROP;                PLPROP=PLPROP*2.7;        OLPROP=OLPROP*3.2;                PLPROP=PLPROP/(PLPROP+OLPROP);        OLPROP=1-PLPROP;               for J=1:12;            DAT(1,J)=(DAT(1,J)-(DAT(2,J)*OLPROP+DAT(3,J)*PLPROP)...                *0.01)/0.99;        end                out = [out;DAT(1,1:12)];                        PLAG_comp = [PLAG_comp; DAT(3,:)];         OL_comp = [OL_comp; DAT(2,:)];                 FractionatingMode = [0 0 OLPROP PLPROP 0 0 0];                       D = FractionatingMode*Dmatrix';        bulkD(:,test+1) = [D];                          DAT_Trace(test+1,:) = (DAT_Trace(test,:) - 0.01.*(DAT_Trace(test,:).*D))./0.99;         FractionatingMode_Out(test+1,:) =FractionatingMode;                                 %fprintf(fid8,'\n');        %fprintf(fid8,'%8.4f',DAT(1,1:12));        if 7.9 < DAT(1,6) & 8.1 > DAT(1,6)            Mg8cnt = Mg8cnt + 1;            Mg8(Mg8cnt,:) = DAT(1,1:12);        end                        [QTZ,PLAG,OLIV,CPX,QQTZ,QOLIV,QCPX,PPLAG,POLIV,PCPX] = Tormey(DAT(1,:));                [HJOLIV,HJPLAG,HJCPX,HJQTZ,QHJQTZ,QHJOLIV,QHJCPX,PHJPLAG,...            PHJOLIV,PHJCPX,MAL2O3,MCAO,MMGO,HJAL2O3,HJCAO,HJMGO] ...            = OPAM (DAT,P);                [XQOPAM,YQOPAM] = PLOT_TERN_TRACE(QHJOLIV,QHJCPX);	           [XMODE,YMODE] = PLOT_TERN_TRACE(QOLIV,QCPX);	            DIFF2=sqrt((XQOPAM-XMODE)^2+(YQOPAM-YMODE)^2);        B=A-DIFF2;        if (B > 0);            A=DIFF2;        end        test = test+1;        OP_first = OP_first+1;         if OP_first ==1             CHANGEINDICIES(1) = test;          end            end%test    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    % OLIVINE, PLAGIOCLASE AND AUGITE FRACTIONATION    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                %SMB: is this because it does an extra step?  should the previous        %step just be deleated? %     for J=1:12;%         DAT(1,J)=(DAT(1,J)+(DAT(2,J)*OLPROP+DAT(3,J)*PLPROP)...%             *0.01)/1.01;%     end%     for J=1:12;%         DAT(1,J)=(DAT(1,J)+(DAT(2,J)*OLPROP+DAT(3,J)*PLPROP)...%             *0.01)/1.01;%     end    %%% 500    [QTZ,PLAG,OLIV,CPX,QQTZ,QOLIV,QCPX,PPLAG,POLIV,PCPX] = Tormey(DAT(1,:));        [HJOLIV,HJPLAG,HJCPX,HJQTZ,QHJQTZ,QHJOLIV,QHJCPX,PHJPLAG,PHJOLIV,...        PHJCPX,MAL2O3,MCAO,MMGO,HJAL2O3,HJCAO,HJMGO] = OPAM_TRACE(DAT,P);        [ROQTZ,ROOLIV,ROPLAG,ROCPX,PROPLAG,PROOLIV,PROCPX,...        QROQTZ,QROOLIV,QROCPX] = OPALM_TRACE(DAT,P);        [QCPXQTZ,QCPXOLIV,QCPXCPX,PCPXPLAG,PCPXOLIV,PCPXCPX,CPXOLIV,CPXPLAG,...        CPXCPX,CPXQTZ] = Tormey(DAT(4,:));    %test = 0;    % FORB = input(['FRACTIONATION = 1, ADDITION = 2 \n']);    FORB = 1;    STEPS = 10;%     %     PROCPX%     PROOLIV%     PROPLAG            OLPROP1 = 0.07;    PLPROP1 = 0.5;    CPXPROP1 = 0.43;            OLPROP2 = 0.07;    PLPROP2 = 0.5;    CPXPROP2 = 0.43;    KDCPX1=DAT(1,6)/DAT(1,5)*DAT(4,5)/DAT(4,6);        %%% 300    %while CONTROL <= 90 || (FORB == 1 & DAT(1,6) > 4.0 & DAT(1,6) < 9.0)    CONTROL=1;    OPC_first = 0;     while CONTROL <= 90 || ( DAT(1,6) > 4.0 && DAT(1,6) < 9.0)        [DAT(2,1),DAT(2,5),DAT(2,6)] = OLCOM_TRACE(DAT,KDOL);	            [DAT(3,1),DAT(3,3),DAT(3,7),DAT(3,8)] = OPAMPL_TRACE(DAT,P);        [DAT(4,1),DAT(4,2),DAT(4,3),DAT(4,5),DAT(4,6),...            DAT(4,7),DAT(4,8),DAT(4,13)] = OPAMCPX_TRACE(DAT,P,KDCPX);        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        %  OUTPUT FRACTIONATION RESULT                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                    if CONTROL < STEPS;            OLPROP=OLPROP1;            PLPROP=PLPROP1;            CPXPROP=CPXPROP1;        else            OLPROP=OLPROP2;            PLPROP=PLPROP2;            CPXPROP=CPXPROP2;        end                                if DAT(1,6) > 3.0                           if FORB == 1;                for J=1:12                    DAT(1,J)=(DAT(1,J)-(DAT(2,J)*OLPROP+DAT(3,J)*PLPROP+...                        DAT(4,J)*CPXPROP)*0.01)/0.99;                end            else                %             for J=1:12                %                 DAT(1,J)=(DAT(1,J)+(DAT(2,J)*OLPROP+DAT(3,J)*PLPROP+...                %                     DAT(4,J)*CPXPROP)*0.01)/1.01;                %             end            end                        FractionatingMode = [CPXPROP 0 OLPROP PLPROP 0 0 0];            D = FractionatingMode*Dmatrix';            bulkD(:,test+1) = [D];            %D = Dmatrix*FractionatingMode';            DAT_Trace(test+1,:) = (DAT_Trace(test,:) - 0.01.*(DAT_Trace(test,:).*D))./0.99;            FractionatingMode_Out(test+1,:) =FractionatingMode;                    else            MagnetiteComp = [0.34	18.5	2.49	0	67.9	0.49	3.9	0	0	0	0	0	0];             %MagnetiteComp = [0.34	21	2.49	0	67.9	0.49	3.9	0	0	0	0	0	0];            MAGPROP = 0.07;            OLPROP=OLPROP*(1-MAGPROP);            PLPROP=PLPROP*(1-MAGPROP);            CPXPROP=CPXPROP*(1-MAGPROP);                        if FORB == 1;                for J=1:12                    DAT(1,J)=(DAT(1,J)-...                        (MagnetiteComp(J)*MAGPROP+...                        DAT(2,J)*OLPROP+...                        DAT(3,J)*PLPROP+...                        DAT(4,J)*CPXPROP)...                        *0.01)/0.99;                end                                FractionatingMode = [CPXPROP 0 OLPROP PLPROP 0 0 MAGPROP];                                  D = FractionatingMode*Dmatrix';                bulkD(:,test+1) = [D];                %D = Dmatrix*FractionatingMode';                DAT_Trace(test+1,:) = (DAT_Trace(test,:) - 0.01.*(DAT_Trace(test,:).*D))./0.99;                FractionatingMode_Out(test+1,:) =FractionatingMode;                            else            end                                end           out = [out;DAT(1,1:12)];                                test=test+1;        OPC_first = OPC_first+1;        if OPC_first ==1            CHANGEINDICIES(2) = test;        end        %fprintf(fid8,' \n');        %fprintf(fid8,'%8.4f',DAT(1,1:12));        if 7.9 < DAT(1,6) & 8.1 > DAT(1,6)            Mg8cnt = Mg8cnt + 1;            Mg8(Mg8cnt,:) = DAT(1,1:12);        end        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        % PLOT RESULT                                                          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        [QTZ,PLAG,OLIV,CPX,QQTZ,QOLIV,QCPX,PPLAG,POLIV,PCPX] = Tormey(DAT(1,:));        [HJOLIV,HJPLAG,HJCPX,HJQTZ,QHJQTZ,QHJOLIV,QHJCPX,PHJPLAG,PHJOLIV,...            PHJCPX,MAL2O3,MCAO,MMGO,HJAL2O3,HJCAO,HJMGO] = OPAM_TRACE(DAT,P);        [ROQTZ,ROOLIV,ROPLAG,ROCPX,PROPLAG,PROOLIV,PROCPX,...            QROQTZ,QROOLIV,QROCPX] = OPALM_TRACE(DAT,P);        [QCPXQTZ,QCPXOLIV,QCPXCPX,PCPXPLAG,PCPXOLIV,PCPXCPX,CPXOLIV,...            CPXPLAG,CPXCPX,CPXQTZ] = Tormey(DAT(4,:));        CONTROL=CONTROL+1;    end    %test    % interpolate to find Mg8 values    % 1    2    3     4     5   6   7   8    9   10   11    12    % SIO2 TIO2 AL2O3 FE2O3 FEO MGO CAO NA2O K2O P2O5 CR2O3 MNO\n']);        Mg8 = out;    %Mg8(isnan(Mg8)) = 0 ;        if any(inputfrommodel) > 0    clear Si8 Ti8 Al8 Fe2O3_8 Fe8 MgO8 Ca8 Na8 K8 P8 Cr8 Mn8 Mg8_int%     Si8 = interp1(Mg8(:,6),Mg8(:,1),8.0);%     Ti8 = interp1(Mg8(:,6),Mg8(:,2),8.0);%     Al8 = interp1(Mg8(:,6),Mg8(:,3),8.0);%     Fe2O3_8 = interp1(Mg8(:,6),Mg8(:,4),8.0);%     Fe8 = interp1(Mg8(:,6),Mg8(:,5),8.0);%     MgO8 = 8.0;%     Ca8 = interp1(Mg8(:,6),Mg8(:,7),8.0);%     Na8 = interp1(Mg8(:,6),Mg8(:,8),8.0);%     K8  = interp1(Mg8(:,6),Mg8(:,9),8.0);%     P8  = interp1(Mg8(:,6),Mg8(:,10),8.0);%     Cr8 = interp1(Mg8(:,6),Mg8(:,11),8.0);%     Mn8 = interp1(Mg8(:,6),Mg8(:,12),8.0);%     Mg8_int = [Si8 Ti8 Al8 Fe2O3_8 Fe8 MgO8 Ca8 Na8 K8 P8 Cr8 Mn8];    % inputfrommodel% [SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,CaO,K2O,Na2O] wt%;    % 1    2    3     4     5   6   7   8    9   10   11    12    % 1SIO2 2TIO2 3AL2O3 4FE2O3 5FEO 6MGO 7CAO 8NA2O 9K2O 10P2O5 CR2O3 MNO\n']);    % %     Mg8_int_out = [Si8 Ti8 Al8 Cr8 Fe8 MgO8 Ca8 K8 Na8  ((MgO8/40.3044))./(((MgO8./40.3044)+((Fe8./71.8464))))];%     out_out = [out(:,1:3) out(:,11) out(:,5:7) out(:,9) out(:,8) ((out(:,6)/40.3044))./(((out(:,6)./40.3044)+((out(:,5)./71.8464))))  ((out(:,6)/40.3044))./(((out(:,6)./40.3044)+((out(:,5)./71.8464))))...%         (out(:,8)+out(:,9))./(out(:,8)+out(:,9)+out(:,7)) out(:,7)./out(:,3)];%         out_out = [out  ((out(:,6)/40.3044))./(((out(:,6)./40.3044)+((out(:,5)./71.8464))))...        (out(:,8)+out(:,9))./(out(:,8)+out(:,9)+out(:,7)) out(:,7)./out(:,3)];    inputfrommodel=[inputfrommodel...        ((inputfrommodel(:,6)/40.3044))./(((inputfrommodel(:,6)./40.3044)+((inputfrommodel(:,5)./71.8464))))...        (inputfrommodel(:,8)+inputfrommodel(:,9))./(inputfrommodel(:,8)+inputfrommodel(:,9)+inputfrommodel(:,7)) ...        inputfrommodel(:,7)./inputfrommodel(:,3)]; if size(inputfrommodel,2)==12;     inputfrommodel(13) = inputfrommodel(:,7)./inputfrommodel(:,3);end  out_out = [   inputfrommodel; out_out];      f= [out_out(:,6)];val = 8; %value to findtmp = abs(f-val);[idx idy] = min(tmp); %index of closest valueclosest = idy;Mg8_int_out(11) = 1-((1-.01)^closest);allindices = [0:size(out_out,1)-1];Fout =1-((1-.01).^(allindices));  PINCall = PINC.*ones(size(Fout)); FractionatingMode_Out = [FractionatingMode_Out PINCall' Fout'];    else                out_out= NaN*ones(1,13);        Mg8_int_out = NaN*ones(1,13);            FractionatingMode_Out = zeros(1,8); DAT_Trace=zeros(1,45);             end%closest = f(idx); %closest value    % %     fprintf(fid6,'\n');% %     fprintf(fid6,['Mg8',file(kk,:),' = [']);% %     fprintf(fid6,'%8.4f',Mg8_int);fprintf(fid6,'];');% %     fprintf(fid6,'\n');% %     clear Mg8% %     fprintf(fid8,'];');% %    fclose(fid8);     %%% 1000 ENDend%fclose(fid6);out_out = [out_out Fout']; % size(out_out)% size(Fout)endendfunction [X,Y] = PLOT_TERN_TRACE (OL,CPX)% Subroutine to be used in conjunction with Yang et al. 1996.% This code has been adapted from FORTRAN by PM Gregg, 2007.%% THIS SUBROUTINE RECALCULATES DATA POINTS ON TERNARIES INTO X-Y COORDINATE% SUBROUTINE PLOT (OL,CPX,X,Y)% REAL OL,CPX,X,YX=(2-CPX/100-2*OL/100)/1.732;Y=CPX/100;endfunction [QTZ,PLAG,OLIV,CPX,QQTZ,QOLIV,QCPX,PPLAG,POLIV,PCPX] = Tormey(DAT)% Subroutine to be used in conjunction with Yang et al. 1996.% This code has been adapted from FORTRAN by PM Gregg, 2007.% It calculates tormey componentsMSIO2=DAT(1)/60.09;MTIO2=DAT(2)/79.9;MAL2O3=DAT(3)/102*2;MFE2O3=DAT(4)/159.7*2;MFEO=DAT(5)/71.85;MMGO=DAT(6)/40.3;MCAO=DAT(7)/56.08;MNA2O=DAT(8)/62*2;MK2O=DAT(9)/94.2*2;MP2O5=DAT(10)/141.9*2;MCR2O3=DAT(11)/152*2;SUM=MSIO2-MCAO-(2*(MNA2O+MK2O))+MTIO2+(2*MP2O5)+MFE2O3;QTZ=(MSIO2-(MMGO+MFEO)/2-(1.5*MCAO)-(MAL2O3/4));QTZ=(QTZ-(2.75*(MNA2O+MK2O)));QTZ=(QTZ+(MCR2O3/2)+(MTIO2/2)+(2.5*MP2O5))/SUM;PLAG=(MAL2O3+MNA2O-MK2O)/(2*SUM);OLIV=((MFEO+MMGO-MTIO2-MCAO)/2+(MAL2O3-MNA2O-MK2O-MCR2O3)...    /4+(0.83333*MP2O5))/SUM;CPX=(MCAO-(MAL2O3/2)+((MK2O+MNA2O)/2)-(1.66667*MP2O5))/SUM;OX=((MCR2O3/2)+MTIO2+MFE2O3)/SUM;AP=MP2O5/(3*SUM);OR=MK2O/SUM;QTZ=QTZ;PLAG=4*PLAG;OLIV=2*OLIV;CPX=3*CPX;QQTZ=100*QTZ/(QTZ+OLIV+CPX);QOLIV=100*OLIV/(QTZ+OLIV+CPX);QCPX=100*CPX/(QTZ+OLIV+CPX);PPLAG=100*PLAG/(PLAG+OLIV+CPX);POLIV=100*OLIV/(PLAG+OLIV+CPX);PCPX=100*CPX/(PLAG+OLIV+CPX);endfunction [ROQTZ,ROOLIV,ROPLAG,ROCPX,PROPLAG,PROOLIV,PROCPX,QROQTZ,...    QROOLIV,QROCPX] = OPALM_TRACE(DAT,P)% Subroutine to be used in conjunction with Yang et al. 1996.% This code has been adapted from FORTRAN by PM Gregg, 2007.%% SUBROUTINE OPALM (DAT,P,ROQTZ,ROOLIV,ROPLAG,ROCPX,PROPLAG,%                   PROOLIV,PROCPX,QROQTZ,QROOLIV,QROCPX)% DIMENSION DAT(51,12)% REAL DAT(51,12),MMGO,MFEO,MG,TIO2,NAK,P% REAL ROQTZ,ROOLIV,ROPLAG,ROCPX% REAL QROQTZ,QROOLIV,QROCPX,PROPLAG,PROOLIV,PROCPXMMGO=DAT(1,6)/40.3;MFEO=DAT(1,5)/71.85;MG=MMGO/(MMGO+MFEO);TIO2=DAT(1,2);NAK=(DAT(1,8)+DAT(1,9))/(DAT(1,8)+DAT(1,9)+DAT(1,7));ROQTZ=0.2-0.012*(P-0.001)+0.00*(1-MG)-0.252*NAK-0.00*TIO2;ROOLIV=0.12+0.008*(P-0.001)+0.229*(1-MG)-0.335*NAK-0.009*TIO2;ROPLAG=0.419+0.012*(P-0.001)-0.135*(1-MG)+0.691*NAK-0.00*TIO2;ROCPX=0.261-0.008*(P-0.001)-0.069*(1-MG)-0.104*NAK+0.012*TIO2;QROQTZ=100*ROQTZ/(ROQTZ+ROOLIV+ROCPX);QROOLIV=100*ROOLIV/(ROQTZ+ROOLIV+ROCPX);QROCPX=100*ROCPX/(ROQTZ+ROOLIV+ROCPX);PROPLAG=100*ROPLAG/(ROPLAG+ROOLIV+ROCPX);PROOLIV=100*ROOLIV/(ROPLAG+ROOLIV+ROCPX);PROCPX=100*ROCPX/(ROPLAG+ROOLIV+ROCPX);endfunction [PLSIO2,PLAL2O3,PLCAO,PLNA2O] = OPAMPL_TRACE(DAT,P)% Subroutine to be used in conjunction with Yang et al. 1996.% This code has been adapted from FORTRAN by PM Gregg, 2007.%% SUBROUTINE OPAMPL_TRACE(DAT,P,PLSIO2,PLAL2O3,PLCAO,PLNA2O)% DIMENSION DAT(51,12)% REAL DAT(51,12),MSIO2,MKALO2,MNAALO2,MCAAL2O4,MCAO,MMGO,MFEO% REAL MTIO2,MPO25,SUM,SUM1,A,MPLAN,MPLAB% REAL MPLCAO,MPLNA2O,MPLAL2O3,MPLSIO2,PLCAO,PLNA2O,PLSIO2,PLAL2O3    MSIO2   = DAT(1,1)/60.09;MKALO2  = DAT(1,9)/94.2*2;MNAALO2 = DAT(1,8)/62*2;MCAAL2O4= (DAT(1,3)/102*2-MKALO2-MNAALO2)/2;MCAO = DAT(1,7)/56-MCAAL2O4;MMGO=DAT(1,6)/40.3;MFEO=DAT(1,5)/71.85;MTIO2=DAT(1,2)/79.9;MPO25=DAT(1,10)/141.9*2;SUM=MSIO2+MKALO2+MNAALO2+MCAAL2O4+MCAO+MMGO+MFEO+MTIO2+MPO25;MSIO2=MSIO2/SUM;MKALO2=MKALO2/SUM;MNAALO2=MNAALO2/SUM;MCAAL2O4=MCAAL2O4/SUM;A=MCAAL2O4/(MNAALO2*MSIO2)*exp(11.1068-0.0338*P-4.4719* ...    (1-MNAALO2)*(1-MNAALO2)-6.9707*(1-MKALO2)*(1-MKALO2));% MPLCAO=(KD*MLIQCAO/MLIQNA2O/2)/((KD*MLIQCAO/MLIQNA2O/2)+1)MPLAN=A/(1+A);MPLAB=1-MPLAN;MPLCAO=MPLAN;MPLNA2O=MPLAB/2;MPLAL2O3=(2*MPLAN+MPLAB)/2;MPLSIO2=2*MPLAN+3*MPLAB;PLCAO=MPLCAO*56.08;PLNA2O=MPLNA2O*62;PLAL2O3=MPLAL2O3*102;PLSIO2=MPLSIO2*60.09;SUM1=PLCAO+PLNA2O+PLSIO2+PLAL2O3;PLCAO=PLCAO/SUM1*100;PLNA2O=PLNA2O/SUM1*100;PLAL2O3=PLAL2O3/SUM1*100;PLSIO2=PLSIO2/SUM1*100;endfunction [HJOLIV,HJPLAG,HJCPX,HJQTZ,QHJQTZ,QHJOLIV,QHJCPX,PHJPLAG,...    PHJOLIV,PHJCPX,MAL2O3,MCAO,MMGO,HJAL2O3,HJCAO,HJMGO] = OPAM_TRACE(DAT,P)% Subroutine to be used in conjunction with Yang et al. 1996.% This code has been adapted from FORTRAN by PM Gregg, 2007.%% SUBROUTINE OPAM (DAT,P,HJOLIV,HJPLAG,HJCPX,HJQTZ,QHJQTZ,QHJOLIV,%                  QHJCPX,PHJPLAG,PHJOLIV,PHJCPX,MAL2O3,MCAO,MMGO,%                  HJAL2O3,HJCAO,HJMGO)% DIMENSION DAT(51,12)% REAL DAT(51,12),P% REAL HJQTZ,HJOLIV,HJPLAG,HJCPX,HJAL2O3,HJCAO,HJMGO% REAL QHJQTZ,QHJOLIV,QHJCPX,PHJPLAG,PHJOLIV,PHJCPX% REAL MSIO2,MTIO2,MAL2O3,MFE2O3,MFEO,MMGO,MCAO,MNA2O% REAL MK2O,MP2O5,MCR2O3,MTOTAL,OX,AP,OR% REAL SUM,QTZ,OLIV,CPX,PLAG,QQTZ,QOLIV,QCPX% REAL PPLAG,POLIV,PCPX% clear P% P = 0.001;MSIO2=DAT(1,1)/60.09;MTIO2=DAT(1,2)/79.9;MAL2O3=DAT(1,3)/102*2;MFE2O3=DAT(1,4)/159.7*2;MFEO=DAT(1,5)/71.85;MMGO=DAT(1,6)/40.3;MCAO=DAT(1,7)/56.08;MNA2O=DAT(1,8)/62*2;MK2O=DAT(1,9)/94.2*2;MP2O5=DAT(1,10)/141.9*2;MCR2O3=DAT(1,11)/152*2;MTOTAL=MSIO2+MTIO2+MAL2O3+MFE2O3+MFEO+MMGO+MCAO+MNA2O+MK2O+MP2O5+MCR2O3;MSIO2=MSIO2/MTOTAL;MTIO2=MTIO2/MTOTAL;MAL2O3=MAL2O3/MTOTAL;MFE2O3=MFE2O3/MTOTAL;MFEO=MFEO/MTOTAL;MMGO=MMGO/MTOTAL;MCAO=MCAO/MTOTAL;MNA2O=MNA2O/MTOTAL;MK2O=MK2O/MTOTAL;MP2O5=MP2O5/MTOTAL;%MCR2O3=CR2O3/MTOTAL;  < -- this line didn't work.  :P changed to next lineMCR2O3=MCR2O3/MTOTAL;HJAL2O3=0.236435+0.002182*P+0.109199*MNA2O+0.593295*MK2O-...    0.349718*MTIO2-0.299419*MFEO-0.130353*MSIO2;HJMGO=-0.276789+0.001139*P-0.542727*MNA2O-0.947294*MK2O-0.117029*MTIO2...	-0.489758*MFEO+2.085584*MSIO2-2.400115*MSIO2^2;HJCAO=1.132807-0.003388*P-0.569004*MNA2O-0.776108*MK2O-0.672285*MTIO2...	-0.213801*MFEO-3.355327*MSIO2+2.830219*MSIO2^2;SUM=MSIO2-HJCAO-(2*(MNA2O+MK2O))+MTIO2+(2*MP2O5)+MFE2O3;QTZ=(MSIO2-(HJMGO+MFEO)/2-(1.5*HJCAO)-(HJAL2O3/4));QTZ=(QTZ-(2.75*(MNA2O+MK2O)));QTZ=(QTZ+(MCR2O3/2)+(MTIO2/2)+(2.5*MP2O5))/SUM;PLAG=(HJAL2O3+MNA2O-MK2O)/(2*SUM);OLIV=((MFEO+HJMGO-MTIO2-HJCAO)/2+(HJAL2O3-MNA2O-MK2O-MCR2O3)/...    4+(0.83333*MP2O5))/SUM;CPX=(HJCAO-(HJAL2O3/2)+((MK2O+MNA2O)/2)-(1.66667*MP2O5))/SUM;OX=((MCR2O3/2)+MTIO2+MFE2O3)/SUM;AP=MP2O5/(3*SUM);OR=MK2O/SUM;HJQTZ=QTZ;HJPLAG=4*PLAG;HJOLIV=2*OLIV;HJCPX=3*CPX;QHJQTZ=100*HJQTZ/(HJQTZ+HJOLIV+HJCPX);QHJOLIV=100*HJOLIV/(HJQTZ+HJOLIV+HJCPX);QHJCPX=100*HJCPX/(HJQTZ+HJOLIV+HJCPX);PHJPLAG=100*HJPLAG/(HJPLAG+HJOLIV+HJCPX);PHJOLIV=100*HJOLIV/(HJPLAG+HJOLIV+HJCPX);PHJCPX=100*HJCPX/(HJPLAG+HJOLIV+HJCPX);endfunction [OLSIO2,OLFEO,OLMGO] = OLCOM_TRACE(DAT,KD)% Subroutine to be used in conjunction with Yang et al. 1996.% This code has been adapted from FORTRAN by PM Gregg, 2007.%% SUBROUTINE OLCOM_TRACE(DAT,KD,OLSIO2,OLFEO,OLMGO)% DIMENSION DAT(51,12)% 	REAL DAT(51,12),KD,MOLMGO,MOLFEO,OLMGO,OLFEO,OLSIO2,SUM% C	WRITE (9,50) (DAT(1,J),J=1,12)% 50      FORMAT (12F7.2)% moles MgO in OlMOLMGO=1/((KD*((DAT(1,5)/71.85)/(DAT(1,6)/40.3)))+1);%MOLMGO=1/((KD*DAT(1,5)/71.85/(DAT(1,6)/40.3))+1)% moles FeO in OlMOLFEO=1-MOLMGO;OLMGO=MOLMGO*40.3;OLFEO=MOLFEO*71.85;OLSIO2=0.5*60.09;SUM=OLMGO+OLFEO+OLSIO2;OLMGO=OLMGO/SUM*100;OLFEO=OLFEO/SUM*100;OLSIO2=OLSIO2/SUM*100;endfunction [CPXSIO2,CPXTIO2,CPXAL2O3,CPXFEO,CPXMGO,CPXCAO,CPXNA2O,TEMP] = ...    OPAMCPX_TRACE(DAT,P,KDCPX)% Subroutine to be used in conjunction with Yang et al. 1996.% This code has been adapted from FORTRAN by PM Gregg, 2007.%% SUBROUTINE OPAMCPX (DAT,P,KDCPX,CPXSIO2,CPXTIO2,CPXAL2O3,%                       CPXFEO,CPXMGO,CPXCAO,CPXNA2O,TEMP)     % DIMENSION DAT(51,12)% REAL DAT(51,12),P,KDCPX,TEMP% REAL MSI,MTI,MAL,MFE2O3,MFE,MMG,MCA,MNA,MK,MP,MCR,MMN% REAL MCPXSI,MCPXTI,MCPXAL,MCPXALT,MCPXALO,MCPXNA% REAL MCPXMGFE,MCPXMG,MCPXFE,MCPXCA% REAL CPXSIO2,CPXTIO2,CPXAL2O3,CPXFEO,CPXMGO,CPXCAO,CPXNA2O% REAL MTOTAL,TOTAL    TEMP=0.4389*DAT(1,7)+4.14*DAT(1,8)-3.103*DAT(1,2);TEMP=TEMP+976;TEMP=TEMP+4.756*P;TEMP=TEMP+21.715*DAT(1,6);TEMP=TEMP+2.509*DAT(1,5);MSI=DAT(1,1)/60.09;MTI=DAT(1,2)/79.9;MAL=DAT(1,3)/102*2;MFE2O3=DAT(1,4)/159.7*2;MFE=DAT(1,5)/71.85;MMG=DAT(1,6)/40.3;MCA=DAT(1,7)/56.08;MNA=DAT(1,8)/62*2;MK=DAT(1,9)/94.2*2;MP=DAT(1,10)/141.9*2;MCR=DAT(1,11)/152*2;MMN=DAT(1,12)/70.94;MTOTAL=MSI+MTI+MAL+MFE2O3+MFE+MMG+MCA+MNA+MK+MP+MCR+MMN;MSI=MSI/MTOTAL*100;MTI=MTI/MTOTAL*100;MAL=MAL/MTOTAL*100;MFE2O3=MFE2O3/MTOTAL*100;MFE=MFE/MTOTAL*100;MMG=MMG/MTOTAL*100;MCA=MCA/MTOTAL*100;MNA=MNA/MTOTAL*100;MK=MK/MTOTAL*100;MP=MP/MTOTAL*100;MCR=MCR/MTOTAL*100;MMN=MMN/MTOTAL*100;MCPXCA=25.044+1.6237*MCA-0.83215*MFE-0.9531*MMG-0.6167*MAL+0.4457*(MNA+MK);% C	MCPXTI=-36.704+0.4661*MCA+0.3654*MFE+0.2590*MMG% C     +         +0.566*MAL+0.3825*(MNA+MK)+0.7905*MTI% C     +         +0.3020*MSI-0.034364*PMCPXTI=-3.6124+0.022523*MCA+0.0339*MFE+0.10776*MAL+0.0426*(MNA+MK)...    +0.3203*MTI-0.02169*P+0.001046*TEMP;% C        MCPXAL=-31.4919+0.5045*MCA+0.3309*MFE+1.6012*MAL% C     +         +1.2832*MTIMCPXNA=-2.2955+0.082*MNA+0.002214*TEMP-0.02024*MMG;MCPXMGFE=-11.051-1.0582*MCA+1.0893*MFE+1.1713*MMG+0.3848*MAL+0.5044*MSI;MCPXSI=116.082-0.31039*MSI-1.374*MTI-1.5822*MAL-0.5334*MFE-0.706*...    MCA-0.4177*MNA+0.2044*P-0.011125*TEMP;% C	MCPXALT=2*MCPXTI% C	MCPXALO=50-MCPXMGFE-MCPXCA-MCPXTI-MCPXNA% C        MCPXALT=50-MCPXSI% C	MCPXAL=MCPXALO+MCPXALTMCPXMG=((KDCPX*MFE/MMG)+1)/MCPXMGFE;MCPXMG=1/MCPXMG;MCPXFE=MCPXMGFE-MCPXMG;MCPXMG=MCPXMGFE*MCPXMG/(MCPXMG+MCPXFE);MCPXFE=MCPXMGFE*MCPXFE/(MCPXMG+MCPXFE);MCPXAL=98.2359-0.96422*MCPXSI-0.7958*MCPXTI-0.997656*MCPXMG-...    1.02176*MCPXFE-1.01482*MCPXCA-0.931973*MCPXNA;% C	MTOTAL=MCPXCA+MCPXTI+MCPXNA+MCPXMGFE+MCPXALO% C	WRITE (9,1001) MTOTAL,MCPXTI,MCPXAL% C1001    FORMAT (3F7.2)% C	MCPXSI=50-MCPXALT-MCPXTICPXSIO2=MCPXSI*60.09;CPXTIO2=MCPXTI*79.9;CPXAL2O3=MCPXAL*102/2;CPXFEO=MCPXFE*71.85;CPXMGO=MCPXMG*40.3;CPXCAO=MCPXCA*56.08;CPXNA2O=MCPXNA*62/2;TOTAL=CPXSIO2+CPXTIO2+CPXAL2O3+CPXFEO+CPXMGO+CPXCAO+CPXNA2O;CPXSIO2=CPXSIO2/TOTAL*100;CPXTIO2=CPXTIO2/TOTAL*100;CPXAL2O3=CPXAL2O3/TOTAL*100;CPXFEO=CPXFEO/TOTAL*100;CPXMGO=CPXMGO/TOTAL*100;CPXCAO=CPXCAO/TOTAL*100;CPXNA2O=CPXNA2O/TOTAL*100;end